# -*- coding: utf-8 -*-
"""
Created on Fri Aug 17 20:23:03 2012

@author: Eduardo
"""


# TODO make this a real package


import numpy as np
import scipy as spy
from spm_ADEM import *


def spm_dx(dfdx,f,t=exp(-10)):
    ''' function [dx] = spm_dx(dfdx,f,t)
% returns dx(t) = (expm(dfdx*t) - I)*inv(dfdx)*f
% FORMAT [dx] = spm_dx(dfdx,f,[t])
% dfdx   = df/dx
% f      = dx/dt
% t      = integration time: (default t = Inf);
%          if t is a cell (i.e., {t}) then t is set to:
%          exp(t - log(diag(-dfdx))
%
% dx     = x(t) - x(0)
%--------------------------------------------------------------------------
% Integration of a dynamic system using local linearization.  This scheme
% accommodates nonlinearities in the state equation by using a functional of
% f(x) = dx/dt.  This uses the equality
%
%             expm([0   0     ]) = (expm(t*dfdx) - I)*inv(dfdx)*f
%                  [t*f t*dfdx]
%
% When t -> Inf this reduces to
%
%              dx(t) = -inv(dfdx)*f
%
% These are the solutions to the gradient ascent ODE
%
%            dx/dt   = k*f = k*dfdx*x =>
%
%            dx(t)   = expm(t*k*dfdx)*x(0)
%                    = expm(t*k*dfdx)*inv(dfdx)*f(0) -
%                      expm(0*k*dfdx)*inv(dfdx)*f(0)
%
% When f = dF/dx (and dfdx = dF/dxdx), dx represents the update from a
% Gauss-Newton ascent on F.  This can be regularised by specifying {t}
% A heavy regularization corresponds to t = -4 and a light
% regularization would be t = 4. This version of spm_dx uses an augmented
% system and the Pade approximation to compute requisite matrix
% exponentials
%
% references:
%
% Friston K, Mattout J, Trujillo-Barreto N, Ashburner J, Penny W. (2007).
% Variational free energy and the Laplace approximation. NeuroImage.
% 34(1):220-34
%
% Ozaki T (1992) A bridge between nonlinear time-series models and
% nonlinear stochastic dynamical systems: A local linearization approach.
% Statistica Sin. 2:113-135.
%
%__________________________________________________________________________
'''
    # defaults
    #----------------------------------------------------------------------

    # if nargin < 3, t = Inf; end

    # t is a regulariser
    #----------------------------------------------------------------------
    #sw = warning('off','MATLAB:log:logOfZero');
    #if iscell(t)
    #    t  = exp(t{:} - log(diag(-dfdx)));

    #warning(sw);

    # use a [pseudo]inverse if all t > TOL
    #======================================================================

    try:
        f.vec
    except:
        f = spm_vec(f)

    dx = spm_vec(f.data)
    
    if min(t) > exp(16):

        dx = -np.linalg.solve(dfdx,f.vec)
        dx.vec[:] = dx

    else:

        # ensure t is a scalar or matrix
        #----------------------------------------------------------------------
        if isvector(t):
            t = np.diag(t)

        # augment Jacobian and take matrix exponential
        #=====================================================================
        # TODO White in C     
        z = np.empty((dfdx.shape[0]+ 1,dfdx.shape[0]+ 1))
        z[0,:] = 0
        for i in xrange(1,dfdx.shape[0]):
            z[i,0] = t * fv[i-1]
            z[i,:] = t * dfdx[i-1,:]            
        dx.vec[:] = np.real(spy.linalg.expm(z))[:,0];
        
    return dx
            
            
        

        
        

class spm_DEM_diff:
    '''function [u dg df] = spm_DEM_diff(M,u)
% evaluates an active model given innovations z{i} and w{i}
% FORMAT [u dgdv dgdx dfdv dfdx] = spm_ADEM_diff(M,u);
%
% M    - generative model
%
% u.v - causal states - updated
% u.x - hidden states - updated
% u.z - innovation (causal state)
% u.w - innovation (hidden states)
% u.a - [active states]
%
% dg.dv, ...  components of the Jacobian in generalised coordinates
%
% The system is evaluated at the prior expectation of the parameters
%__________________________________________________________________________
'''

    def __init__(self,M):
        
        if M.is_ADEM:
            self.is_ADEM = True
        else:
            self.is_ADEM = False
            u.a      = u.v;
            M.a      = [np.zeros((0,1))]*ml
            M.k      = [0]*ml;
            u.a      = [np.zeros((0,1))]*ml;

        # number of states and parameters
        #======================================================================
        ml    = M.ml               # number of levels
        nv    = sum(M.l)           # number of v (causal states)
        na    = sum(M.k)           # number of a (active states)
        nx    = sum(M.n)           # number of x (hidden states)
 
        # order parameters (n = 1 for static models)
        #----------------------------------------------------------------------
        n     = M.E.n + 1          # order of embedding

 
        # initialise arrays for hierarchical form
        #-----------------------------------------------------------------------
        
        dfdvi = [[]]*ml
        dfdxi = [[]]*ml
        dfdai = [[]]*ml
        dgdvi = [[]]*ml
        dgdxi = [[]]*ml
        dgdai = [[]]*ml
        
        for i in xrange(ml):
            dgdvi[i] = [np.zeros(M.l[i],M.l[i])]
            dgdxi[i] = [np.zeros(M.l[i],M.n[i])]
            dgdai[i] = [np.zeros(M.l[i],M.k[i])]
            dfdvi[i] = [np.zeros(M.l[i],M.l[i])]
            dfdxi[i] = [np.zeros(M.n[i],M.n[i])]
            dfdai[i] = [np.zeros(M.l[i],M.k[i])]

 
        # partition states {x,v,z,w} into distinct vector arrays v{i}, ...
        #----------------------------------------------------------------------
        
        uvv = spm.vec(u.v[0])
        uxv = spm.vec(u.x[0])
        uav = spm.vec(u.a[0])
        uzv = spm.vec(u.z[0])
        uwv = spm.vec(u.w[0])
        
        Mvv = spm.vec(M.v)
        Mxv = spm.vec(M.x)
        Mav = spm.vec(M.a)
        Mzv = spm.vec(M.v)
        Mxv = spm.vec(M.x)
        
        uvv.vec[:] = Mvv.vec
        uxv.vec[:] = Mxv.vec
        uav.vec[:] = Mav.vec
        uzv.vec[:] = Mzv.vec
        uwv.vec[:] = Mxv.vec
        
        self.vi = uvv.data
        self.xi = uxv.data
        self.ai = uave.data
        self.zi = uzv.data
        self.wi = uwv.data
        class d:
            def __init__(self):
                self.da = None
                self.dv = None
                self.dx = None
        self.dg = d()
        self.df = d()

        
    
    def next(self,u):
        
        # Derivatives for Jacobian
        #======================================================================
        # TODO what's that line?        
        #vi{nl} = zi{nl};
        for  i in xrange(ml-1,-1,-1):
               # evaluate
               #---------------------------------------------------------------
               if self.is_ADEM:
                   dgdx,g   = spm_diff(f=M.g[i],n=[0],x=[xi[i],v[i+1],ai[i+1]])
                   dfdx,f   = spm_diff(f=M.f[i],n=[0],x=[xi[i],v[i+1],ai[i+1]])
                   dgdv     = spm_diff(f=M.g[i],n=[1],x=[xi[i],v[i+1],ai[i+1]])[0]
                   dfdv     = spm_diff(f=M.f[i],n=[1],x=[xi[i],v[i+1],ai[i+1]])[0]
                   dgda     = spm_diff(f=M.g[i],n=[2],x=[xi[i],v[i+1],ai[i+1]])[0]
                   dfda     = spm_diff(f=M.f[i],n=[2],x=[xi[i],v[i+1],ai[i+1]])[0]
               else:
                   dgdx,g   = spm_diff(f=M.g[i],n=[0],x=[xi[i],v[i+1],ai[i+1]])
                   dfdx,f   = spm_diff(f=M.f[i],n=[0],x=[xi[i],v[i+1],ai[i+1]])
                   dgdv     = spm_diff(f=M.g[i],n=[1],x=[xi[i],v[i+1],ai[i+1]])[0]
                   dfdv     = spm_diff(f=M.f[i],n=[1],x=[xi[i],v[i+1],ai[i+1]])[0]
                   dgda     = []
                   dfda     = []
                   
               # g(x,v) & f(x,v)
               #---------------------------------------------------------------
               gi[i]    = g
               fi[i]    = f
               vi[i]    = gi[i] + zi[i]
            
               # and partial derivatives. If necesary construct the block diagonal
               # matrices
               #---------------------------------------------------------------
               if isinstance(self.dgdxi,list):
                   self.dgdxi[i][:]  = dgdx
                   self.dgdvi[i]     = np.concatenate((dgdvi,dgdv),axis=1)
                   self.dgdai[i]     = np.concatenate((dgdai,dgda),axis=1)
                   self.dfdxi[i][:]  = dfdx
                   self.dfdvi[i]     = np.concatenate((dfdvi,dfdv),axis=1)
                   self.dfdai[i]     = np.concatenate((dfdai,dfda),axis=1)
               else:
                   # fill the contenst of a block diagonal matrix
                   #-----------------------------------------------------------
                   self.dgdxi[i][:]  = dgdx
                   self.dgdvi[i][:]  = np.concatenate((dgdvi,dgdv),axis=1)
                   self.dgdai[i][:]  = np.concatenate((dgdai,dgda),axis=1)
                   self.dfdxi[i][:]  = dfdx
                   self.dfdvi[i][:]  = np.concatenate((dfdvi,dfdv),axis=1)
                   self.dfdai[i][:]  = np.concatenate((dfdai,dfda),axis=1)
                   
        # concatenate hierarchical arrays
        #--------------------------------------------------------------------------
        self.dg.da[:] = np.linalg.block_diag(self.dgdai)
        self.dg.dv[:] = np.linalg.block_diag(self.dgdvi)
        self.dg.dx[:] = np.linalg.block_diag(self.dgdxi)
        self.df.da[:] = np.linalg.block_diag(self.dvdai)
        self.df.dv[:] = np.linalg.block_diag(self.dvdvi)
        self.df.dx[:] = np.linalg.block_diag(self.dvdxi)
 
        # update generalised coordinates
        #--------------------------------------------------------------------------
        u.v[0]  = vi.vec
        u.x[1]  = fi.vec + u.w[i]
        for i in xrange(1,n-1):
            u.v[i]     = dg.dv*u.v[i].vec + dg.dx*u.x[i].vec + dg.da*u.a[i].vec + u.z[i].vec
            u.x[i + 1] = df.dv*u.v[i].vec + df.dx*u.x[i].vec + df.da*u.a[i].vec + u.w[i].vec
        
        return u



def spm_DEM_embed(self,Y,n,t,dt=1,d=0):
    '''function [y] = spm_DEM_embed(Y,n,t,dt,d)
% temporal embedding into derivatives
% FORMAT [y] = spm_DEM_embed(Y,n,t,dt,d)
%__________________________________________________________________________
% Y    - (v x N) matrix of v time-series of length N
% n    - order of temporal embedding
% t    - time  {bins} at which to evaluate derivatives (starting at t = 1)
% dt   - sampling interval {secs} [default = 1]
% d    - delay (bins) for each row of Y
%
% y    - {n,1}(v x 1) temporal derivatives   y[:] <- E*Y(t)
%==========================================================================
'''

    # get dimensions
    #----------------------------------------------------------------------
    q,N    = Y.shape
    y      = [np.zeros((q,1))] * n
    
    #----------------------------------------------------------------------
    if q == 0:
        return

    # loop over channels
    #----------------------------------------------------------------------
    try:
        iter(d)
    except:
        d = [d]

    def nott(x):
       if x == 0: 
           return 1
       else: 
           return 0
        
    for p in xrange(d):

       # boundary conditions
       #-------------------------------------------------------------------
       s      = (t - d[p])/dt
       k      = np.arange(n)  + np.fix(s - (n + 1.)/2.)
       x      = s - min(k) + 1.
       i      = k < 1
       k      = k*np.logical_not(i) + i
       i      = k > N
       k      = k*np.logical_not(i) + i*N;

       # TODO rewrite in c?

       # Inverse embedding operator (T): cf, Taylor expansion Y(t) <- T*y[:]
       #------------------------------------------------------------------
        
       T = np.zeros((n,n))
        
       for i in xrange(n):
           for j in xrange(n):
               T[i,j] = ((i - x)*dt)^(j - 1)/np.prod(np.arange(j-1)+1);


       # embedding operator: y[:] <- E*Y(t)
       #------------------------------------------------------------------
       E     = np.inv(T)

       E     = np.asmatrix(E)
       Y     = np.asmatrix(Y)
       

       # embed
       #------------------------------------------------------------------
       
       if len(d) == q:
           for i in xrange(n):
               y[i][p,:] = Y[p,k]*E[i,:].T
       else:
           for i in xrange(n):
               y[i]      = Y[:,k]*E[i,:].T

       return y




class spm_DEM_int:
    ''' function [V,X,Z,W] = spm_DEM_int(M,z,w,c)
% Integrates/evaluates a hierarchical model given innovations z{i} and w{i}
% FORMAT [V,X,Z,W] = spm_DEM_int(M,z,w,c);
%
% M{i}    - model structure
% z{i}    - innovations (causes)
% w{i}    - innovations (states)
% c{i}    - exogenous causes
%
% V{i}    - causal states (V{1} = y = response)
% X{i}    - hidden states
% Z{i}    - fluctuations (causes)
% W{i}    - fluctuations (states)

'''
    def __init__(self,M,z,w,c):
        
        # set model indices and missing fields
        #----------------------------------------------------------------------
    
        M.check()

        # innovations
        #----------------------------------------------------------------------
    
        #try, z = spm_cat(z(:)); end
        #try, w = spm_cat(w(:)); end
        #try, c = spm_cat(c(:)); end

        # number of states and parameters
        #----------------------------------------------------------------------
        nt   = z.shape[1]                       # number of time steps
        ml   = M.ml                             # number of levels
        nv   = sum(M.l)                         # number of v (casual states)
        nx   = sum(M.n)                         # number of x (hidden states)

        # order parameters (n= 1 for static models)
        #======================================================================
        dt   = M.E.dt                           # time step
        n    = M.E.n + 1                        # order of embedding
        nD   = M.E.nD                           # number of iterations per sample
        td   = dt/nD                            # integration time for D-Step

        # initialize cell arrays for derivatives z{i} = (d/dt)^i[z], ...
        #----------------------------------------------------------------------
        
        class u_c:
            def __init__(self,n,nv,nx):
                self.v = [np.zeros((nv,1))] * n
                self.x = [np.zeros((nx,1))] * n
                self.z = [np.zeros((nv,1))] * n
                self.w = [np.zeros((nx,1))] * n
        
        u = u_c(n,nv,nx)

        # hyperparameters
        #----------------------------------------------------------------------

        class ph_c:
            def __init__(self,M):
                self.h = M.hE
                self.f = M-gE
        
        ph = ph_c(M)
        
        # initialize with starting conditions
        #----------------------------------------------------------------------
        vi     = M.v
        xi     = M.x
        
        u.v[0] = spm_vec(vi);
        u.x[0] = spm_vec(xi);


        # derivatives for Jacobian of D-step
        #----------------------------------------------------------------------
        Dx    = np.kron(np.eye(n,k=1),np.eye(nx,0))
        Dv    = np.kron(np.eye(n,k=1),np.eye(nv,0))
        # D     = spm_cat(spm_diag({Dv,Dx,Dv,Dx}));
        D     = spy.linalg.block_diag(Dv,Dx,Dv,Dx)
        dfdw  = np.kron(np.eye(n,n),np.eye(nx,nx))

        # initialize conditional estimators of states to be saved (V and X)
        #----------------------------------------------------------------------
        for i in xrange(ml):
            V[i] = np.zeros((M.l[i],nt))
            X[i] = np.zeros((M.n[i],nt))
            Z[i] = np.zeros((M.l[i],nt))
            W[i] = np.zeros((M.n[i],nt))


        # method for state-dependent precision
        #----------------------------------------------------------------------
        
        phv = spm.vec(M.ph)
        pgv = spm.vec(M.pg)
        
        if sum(phv.vec.shape) == 0 and np.sum(pgv.vec.shape):
            state_dependent = 0
            Sz = 1
            Sw = 1
        else:
            state_dependent = 1

        # iterate over sequence (t) and within for static models    
        #======================================================================
        for t in xrange(nt):
            for iD in xrange(nD):
                
                # Get generalised motion of random fluctuations
                #==============================================================

                # sampling time
                #--------------------------------------------------------------
                
                ts   = (t + (iD - 1)/nD)*dt

                # evaluate state-dependent precision
                #--------------------------------------------------------------
                if state_dependent:
                    
                    pu.x = [spm.vec(xi[1:-2])]
                    pu.v = [spm.vec(xi[2:-1])]
                    # TODO implement this
                    p    = spm_LAP_eval(M,pu,ph)
                    Sz   = np.zeros(np.diag(np.exp(-p.h/2))) 
                    Sw   = np.zeros(np.diag(np.exp(-p.g/2))) 


                # derivatives of innovations (and exogenous input)
                #--------------------------------------------------------------
                u.z  = spm_DEM_embed(Sz*z + c,n,ts,dt)
                u.w  = spm_DEM_embed(Sw*w,    n,ts,dt)


                # Evaluate and update states
                #==============================================================

                # evaluate functions
                #--------------------------------------------------------------
                u,dg,df = spm_DEM_diff(M,u);

        # tensor products for Jacobian
        #------------------------------------------------------------------
        dgdv = np.kron(np.eye(n,1),dg.dv)
        dgdx = np.kron(np.eye(n,1),dg.dx)
        dfdv = spy.lin_alg.block_diagonal([[df.dv]]*n)
        dfdx = spy.lin_alg.block_diagonal([[df.dx]]*n)

        # Save realization
        #==================================================================
        vi   = u.v[0]
        xi   = u.x[0]
        zi   = u.z[0]
        wi   = u.w[0]
        
        if iD == 1:
            for i in xrange(nl):
                if M.l[i]: 
                    V[i][:,t] = spm.vec(vi.data[i]).vec
                if M.n[i]: 
                    X[i][:,t] = spm.vec(xi.data[i]).vec
                if M.l[i]: 
                    Z[i][:,t] = spm.vec(zi.data[i]).vec
                if M.n[i]: W[i][:,t] = spm.vec(wi.data[i]).vec

        # Jacobian for update
        #------------------------------------------------------------------
        J = self.jacobian()

        # update states u = {x,v,z,w}
        #------------------------------------------------------------------
        du     = spm_dx(J,D*spm_vec(u),td);

        # and unpack
        #------------------------------------------------------------------
        u      += du
        
    def jacobian(self):
        # TODO
        '''
        J      = spm_cat([dgdv dgdx Dv  []  ;
                          dfdv dfdx []  dfdw;
                          []   []   Dv  []  ;
                          []   []   []  Dx]);
        '''
        return []
        


    
    


def spm_sqrtm(V):
    '''Matrix square root for sparse symmetric positive semi-definite matrices
% FORMAT [K] = spm_sqrtm(V)
%
% This routine covers and extends sqrtm functionality by using a
% computationally expedient approximation that can handle sparse
% symmetric positive semi-definite matrices.
    '''
    # TODO this is actually stupid but is there to keep the interface. Maybe 
    # one day there will be stuff for sparse matrices
    
    if not V.shape[0] == V.shape[1]:
        raise TypeError(msg='Non square matrix')
        
    U,S = np.svd(V)
    np.sqrt(S,S)
    return np.dot(np.dot(U,S),U.T)

class spm_DEM:
    
    def __init__(self,M):
        
        if not isinstance(M,'Model'):
            raise TypeError(msg='M is not a model')
        
        M.check()        
        self.M = M

        self.Y = None

        class p:
            '''Posteriors'''
            def __init__ (self):
                class U:
                    ''' Time depending states'''
                    def __init__(self):
                        self.v = None
                        self.x = None
                        self.e = None
                self.U = U()
                self.P = None
                self.H = None
                
        self.p = p()




        
class spm_DEM_z:
    
    def __init__(self,M,N):
        ''' function [z,w] = spm_DEM_z(M,N)
% creates hierarchical innovations for generating data
% FORMAT [z w] = spm_DEM_z(M,N)
% M    - model structure
% N    - length of data sequence
%
% z{i} - innovations for level i (N.B. z{end} corresponds to causes)
% w{i} - innovations for level i (state noise)
%
% If there is no fixed or hyper parameterized precision, then unit noise is
% created. It is assumed that this will be later modulated by state
% dependent terms, specified by M.ph and M.pg in spm_DEM_int
%__________________________________________________________________________
'''
         
        # temporal convolution matrix (with unit variance)
        #---------------------------------------------------------------------
         
        M.check()
        
        s  = M.E.s + exp(-16)
        dt = M.E.dt 
        t  = np.range(N,dtype='d') * dt
        t  = np.exp(np.power(-f,2)/(2*s^2))
        K  = np.linalg.toeplitz(t)
        K  = np.dot(np.diag(np.power(np.diagonal(np.dot(K,K.T),-0.5)),K))
        
        # create innovations z{i} and w{i}
        #----------------------------------------------------------------------

        z = [] * M.ml
        w = [] * M.ml
        
        for i in xrange(M.ml):
            # causes: assume i.i.d. if precision (P) is zero
            #------------------------------------------------------------------
            P = M.V[i]
            for j in xrange(len(M.Q[i])):
                P = P + np.dot(M.Q[i][j],np.exp(M.hE[i][j]))
            if not np.norm(P,1):
                P = 1
            
            # z{i}  = spm_sqrtm(inv(P))*randn(M(i).l,N)*K;
            
            U,S  = np.svd(P)
            tS   = np.diag(np.power(np.diag(S),-0.5))
            isP  = np.dot(np.dot(U.T,tS),U)
            
            z[i] = np.dot(np.dot(isP,spy.random.randn(M.l[i],N)),K)
            
            # states
            #------------------------------------------------------------------
            
            P = M.W[i]
            for j in xrange(len(M.R[i])):
                P = P + M.R[i][j] * exp(M.gE[i][j])                
            
            if not P == []:
                if np.norm(P,1) == 0:
                    P = 1
                U,S  = np.svd(P)
                tS   = np.diag(np.power(np.diag(S),-0.5))
                isP  = np.dot(np.dot(U.T,tS),U)
                
                #w{i} = spm_sqrtm(inv(P))*randn(M(i).n,N)*K*dt;
                
                w[i] = np.dot(np.dot(isP,spy.random.randn(M.n[i],N)),K) * dt
            else:
                w[i] = np.zeros((0,0))
            
        self.z = z
        self.w = w
        
            
class spm_DEM_generate:
    ''' function [DEM] = spm_DEM_generate(M,U,P,h,g)
    % Generates data for a Hierarchical Dynamic Model (HDM)
    % FORMAT [DEM] = spm_DEM_generate(M,N,P,h,g): N-samples using z
    % FORMAT [DEM] = spm_DEM_generate(M,U,P,h,g): size(U,2) samples using U
    %
    % M(i)     - HDM
    % U(n x N} - causes or N number of causes
    % P{i}     - model-parameters for level i (defaults to M.pE)
    % h{i}     - log-precisions   for level i (defaults to 32 - no noise)
    % g{i}     - log-precisions   for level i (defaults to 32 - no noise)
    %
    % generates
    % DEM.M    - hierarchical model (checked)
    % DEM.Y    - responses or data
    %
    % and true causes NB: v{end} = U or z{end} (last level innovations)
    % DEM.pU.v 
    % DEM.pU.x
    % DEM.pU.e
    % DEM.pP.P
    % DEM.pH.h
    '''
    
    def __init__(self,M,U,P=None,h=None,g=None):
        
        if not isinstance(M,'Model'):
            raise TypeError(msg='M is not a model')

        # check model
        #----------------------------------------------------------------------
        M.check()
        self.DEM = spm_DEM(M)
        
        # sequence length specified by priors on causes
        #----------------------------------------------------------------------
        
        try:
            U.shape
        except AttributeError:
            raise TypeError(msg='please check U')
            
        if np.prod(U.shape) > 1:
            N = U.shape[1]
        else:
            N = U
            U = np.zeros(M.l[-1],N)
            
 
        # initialize model-parameters if specified
        #----------------------------------------------------------------------
         
        if P is None:
            P = [M.pE]
        elif not isinstance(P,list):
            P = [P]
        if h is None:
            h = []
        elif not isinstance(h,list):
            h = [h]
        if g is None:
            g = []
        elif not isinstance(h,list):
            g = [h]
 
        # transcribe parameters and hyperparameters into prior expectations
        #----------------------------------------------------------------------
        ml     = M.ml
        
        for i in xrange(ml):
            try:
                tP  = spm.vec(P[i])
                tpE = spm.vec(M.pE[i])
                if len(spm.vec(P[i]).vec) == len(spm.vec(M.pE[i])):
                    tpE.vec[:] = tP.vec
                    M.pE[i]    = tpE.data
            except:
                try:
                    M.pE[i] = P[i]
                except:
                    pass
        for i in xrange(ml):
            try:
                M.hE[i] = h[i]
            except:
                
                M.hE[i] = (M.hE[i] - M.hE[i]) + 32
                
        for i in xrange(ml):
            try:
                M.gE[i] = g[i]
            except:
                M.gE[i] = (M.gE[i] - M.gE[i]) + 32
            
        
        # re-set M and create innovations
        #----------------------------------------------------------------------

        M.check()        
        
        rn = spm_DEM_z(M,N)
        
        
        # place exogenous causes in cell array
        #----------------------------------------------------------------------
        u = []*ml
        for i in xrange(ml):
            u[i] = np.zeros(M.l[i],N)
        
        u[ml] = U

        # integrate HDM to obtain causal (v) and hidden states (x)
        #----------------------------------------------------------------------
        # [v,x,z,w] = spm_DEM_int(M,z,w,u);
        
        DEM_int = spm_DEM_int(M,rn,u)
        
        
        # Fill in DEM with response and its causes
        #----------------------------------------------------------------------
        self.DEM.Y      = v[0]
        self.DEM.p.U.v  = v;
        self.DEM.p.U.x  = x;
        self.DEM.p.U.z  = z;
        self.DEM.p.U.w  = w;
        self.DEM.p.P.P  = {M.pE};
        self.DEM.p.H.h  = {M.hE};
        self.DEM.p.H.g  = {M.gE};
